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Abstract —Multi-target tracking is an important problem in 
civilian and military applications. This paper investigates multi¬ 
target tracking in distributed sensor networks. Data association, 
which arises particularly in multi-object scenarios, can be tack¬ 
led by various solutions. We consider sequential Monte Carlo 
implementations of the Probability Hypothesis Density (PHD) 
filter based on random finite sets. This approach circumvents 
the data association issue by jointly estimating all targets in 
the region of interest. To this end, we develop the Diffusion 
Particle PHD Filter (D-PPHDF) as well as a centralized version, 
called the Multi-Sensor Particle PHD Filter (MS-PPHDF). Their 
performance is evaluated in terms of the Optimal Subpattern 
Assignment (OSPA) metric, benchmarked against a distributed 
extension of the Posterior Cramer-Rao Lower Bound (PCRLB), 
and compared to the performance of an existing distributed PHD 
Particle Filter. 

Index Terms —Multi-target tracking, distributed target track¬ 
ing, Particle Filter, PHD Filter, autonomous sensor networks, 
Posterior Cramer-Rao Lower Bound 


I. Introduction 

T HE problem of multi-target tracking (MTT) is becoming 
increasingly important in many military and civilian 
applications such as air and ground traffic control, harbor 
surveillance, maritime traffic control, or video communication 
and surveillance 0-0. Distributed sensor networks offer 
a desirable platform for MTT applications due to their low 
cost and ease of deployment, their lack of a single point 
of failure, as well as their inherent redundancy and fault- 
tolerance 0 . A comprehensive overview of the state-of-the- 
art of distributed single-target tracking (STT) is given in 0. 
Distributed versions of the Kalman Filter 0, 0 and its 
nonlinear, non-Gaussian counterpart, the Particle Filter (PF) 
0, have been well-studied. However, they cannot be applied 
directly to MTT as they do not account for the problem of 
data association. Although there are methods such as the 
Joint Probabilistic Data Association Filter (JPDAF) 0 or 
the Multiple Hypothesis Tracker (MHT) 0 that address this 
problem in STT algorithms, the resource constraints in sensor 
networks might pose a challenge on finding suitable distributed 
implementations m. The Probability Hypothesis Density 
(PHD) filter ifTTj . fl2l . in contrast, resorts to the concept of 
random finite sets (RFSs) to circumvent the problem of data 
association altogether. 

In this work, we investigate distributed MTT in a sensor 
network with 1-coverage of the region of interest (ROI), i.e., 
the sensor nodes have non- or barely overlapping fields of view 
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(FOVs) and are distributed such that maximum area coverage 
is attained CCD. An exemplary network layout with these 
properties is depicted in Figure [Tall. Autonomous distribution 
algorithms for realizing such a topology have been studied in 
our previous work 03 - The nodes in the network commu¬ 
nicate with their neighbors in order to collaboratively detect 
and track targets in the ROI. In addition, all of the sensors 
are equipped with a signal processing unit, allowing them to 
form decisions without a fusion center. That way, the network 
can autonomously react to events such as the detection of an 
intruder without relying on a network operator. For the sake of 
simplicity, the network is considered to be static. However, the 
consideration of mobile sensor nodes would enable reactions 
such as target pursuit or escape. 

Since the FOV and communication radius of each node are 
limited, a target is only seen by a subset of the network, which 
changes as the target moves through the ROI. Hence, at each 
time instant, there is an active and an inactive part of the 
network. The goal, thus, is to detect and observe the target 
in a distributed and collaborative fashion as it travels across 
the ROI, rather than reaching a network-wide consensus on its 
state and have the estimate available at each node. 

In the sequel, we develop a distributed Particle PHD filter 
called Diffusion Particle PHD Filter (D-PPHDF), which uses 
neighborhood communication to collaboratively estimate and 
track a single-sensor PHD at each node in the active subnet¬ 
work. In addition, we formulate the Multi-Sensor Particle PHD 
Filter (MS-PPHDF), a centralized extension of the D-PPHDF. 
The performance of both algorithms is evaluated in terms 
of the Optimal Subpattern Assignment (OSPA) metric na, 
which is calculated for the joint set of target state estimates of 
the active subnetwork. Furthermore, a distributed version of 
the Posterior Cramer-Rao Lower Bound (PCRLB) lfT6fl - lfl8ll 
- again averaged over the active subnetwork - is introduced 
and used as a benchmark. 

Other distributed solutions for MTT in a multi-sensor setup 
using the PHD filter have been studied, e.g., in CGD, m, 
El. Contrary to our approach, they either assume overlapping 
FOVs or employ a pairwise communication scheme. The 
common idea, however, is to extend the single-sensor PHD 
filter to the multi-sensor case through communication between 
multiple nodes, or nodes and a fusion center. A more rigorous 
approach for MTT with multiple sensors is to use a multi¬ 
sensor PHD filter l22l , l23l . which seeks to estimate and 
track a single multi-sensor PHD instead of multiple single¬ 
sensor PHDs. In this work, we compare our methods to the 
approach in |20j (adapted to our scenario), which is also based 
on single-sensor PHDs. The consideration of methods based 
on a multi-sensor PHD will be the focus of future work. 
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Fig. 1: (a) Distributed sensor network with 1-coverage of the region of interest and 3 exemplary target tracks, (b) Example 
of tracking 3 targets with the Diffusion Particle PHD Filter (D-PPHDF). The small colored dots represent the target location 
estimates obtained by the respective node with the same color. The light gray dots show the collective measurements obtained 
by all nodes in the network. 


The paper is organized as follows: Section 2 presents the 
considered state-space model and recapitulates the theory of 
RFSs as well as the PHD and the PHD filter. The problem 
of distributed MTT is addressed in Section 3. Here, we 
will first detail our modification of Adaptive Target Birth 
(ATB) before formulating the D-PPHDF and investigate its 
computational complexity and communication load. In Section 
4, the MS-PPHDF is developed and analyzed in terms of 
computational complexity and communication load. Section 
5 is dedicated to simulations. First, the Distributed Posterior 
Cramer-Rao Lower Bound (DPCRLB) is introduced. Then, we 
present the simulation setup and discuss our results. Finally, a 
conclusion is given in Section 6. 

II. Models and Theory 
A. State-Space Model 

A linear state-space model is considered for each target 
at time instant * > 0. The target state vector s tgt (i) = 
[x tgt (i), i: tst (*)] T contains the target location vector tc tgt as 
well as the velocity vector x tgt . For the sake of simplicity, 
we restrict ourselves to a 2D-environment. The target state 
evolves according to the state equation (24): 

s tgt («) = — 1 ) + G(*)n tgt (i). ( 1 ) 

The matrices F and G as well as the vector n tgt will be 
explained shortly. Node k obtains a measurement Zk of the 
target location as given by the measurement equation lf24l : 

z k (i) = H k {i)s t& {i) + keM ( 2 ) 

with M = {m G {1,...,N} \ \\x m (i) - tc tgt (i )|| 2 < R sen } 
denoting the set of all nodes m that are located such that the 


Euclidean distance \\x m (i) — £c tgt (z) ||2 between their location 
x m and the target location x tgt is not greater than their 
sensing radius R sen . Note that N is the total number of nodes 
in the network. Furthermore, n tgt («) ~ A/"(02 ,i, Q(i)) and 
~ M( 0 2 tl ,R k (i)) denote the state and measurement 
noise processes, respectively, with the zero-mean vector O 2.1 = 
[0,0] T . Both noise processes are spatially and temporally 
white, as well as uncorrelated with the initial target state s tgt (0) 
and each other for all i. For the sake of simplicity, we choose 
a time-invariant measurement noise covariance matrix 


Rk(i) = Rk = 0II2, (3) 


where cr^ is the variance of each component of the measure¬ 
ment noise and I n denotes the identity matrix of size n. 

In target tracking, the model matrices are usually chosen to 
be time-invariant and given by ED 


F = 


I2 

£•2,2 


AH2 

I2 


G = 


A £ 

Ail 2 


At r 

— J 2 


Q 


&ql 2 , 


(4) 

where 02,2 is the 2x2 zero matrix. Furthermore, A* is 
the time step interval in seconds with which the state-space 
model progresses. In addition, a% denotes the variance of 
a state noise component. We assume that the sensor nodes 
only obtain location information, for instance, through distance 
and bearing measurements from which an estimate of the 
target location can easily be calculated. This yields a general 
measurement matrix Hk of the form 


H k = [I2 02,2] ■ 


(5) 


B. Random Finite Sets (RFSs) 

A random finite set (RFS) is an unordered finite set that 
is random in the number of its elements as well as in 
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their values (25). Therefore, RFSs are a natural choice for 
representing the multi-target states and measurements in MTT: 
the state and measurement vectors of all targets are collected in 
corresponding RFSs (26), E 7 ). Given the realization —,- i of 
the RFS Ej_i at time instant i — 1, the multi-target state of our 
tracking problem can be described by the RFS E, according 
to 

Ej = U Bi, (6) 

where the survival set denotes the RFS of targets 

that already existed at time step i — 1 and have not exited 
the ROI, i.e., the region covered by the sensor network, in the 
transition to time step i. In addition, the birth set 13, is the 
RFS of new targets that spontaneously appear at the border of 
the ROI at time instant i ID, 02, El- Note that the statistical 
behavior of E, can be described by the conditional probability 
fi\i—l ( l —I 1 —'i—1 )• 

The multi-target measurement model is given by the RFS 

S; as 

Sj = ®i(Sj) U (7) 

where ©i(S,) is the RFS of measurements generated by E,. 

In addition, the RFS Ci (Sj) represents clutter or false alarms. 
Given a realization E, of XI,, the statistical behavior of the 
RFS Si is described by the conditional probability /,(E,;|E,;). 

C. The Probability Hypothesis Density (PHD) 

In analogy to the single-target case, the optimal Bayesian 
filter for MTT recursively propagates the multi-target posterior 
/i|i(5i|S 0 ;i) over time, according to 

, iy \ /i(S i |a i )/ i | i _i(a i |S 0:i —i) 

° :ij “ // i (S i |S)/ i | i _ 1 (S|Eo :i _ 1 K(dS) ( 

fi\i— 1 C 1 —'i |Eo:i— l) — ^ /i| j —1 l|i—1 (^1 So:i —l)/X s (d^.) 

(9) 

where p s is a dominating measure as described in G3- 
This approach requires the evaluation of multiple integrals, 
which makes it even more computationally challenging than 
its single-target counterpart. A common solution is to find a set 
of statistics, e.g., the moments of first or second order, which 
yield a good approximation of the posterior, and propagate 
them instead m. 

The Probability Hypothesis Density (PHD) A|i( s (*)l^ 0 :j) 
is an indirect first-order moment of /ju(Sj|Eo : j) |28j|. It is 
given by the following integral H2, G2: 

A|i(s(*)|£ 0 :i) 

= E [ S ( S ( i ) ~ S n\i))fi[d S ^)\^0:i)ds{i), (10) 

where f f(Y)8Y denotes a set integral. 

The PHD has the following two properties (28) : 

1) The expected number of targets N tgt (i) at time step i is 
obtained by integrating the PHD according to 

Ntgt(i) = J D x \i(s(i)\Y, 0 ,i)ds(i). (11) 


2 ) 


This is in contrast to probability density functions 
(PDFs), which always integrate to 1. 

Estimates of the individual target states can be found 


by searching for the 


N. 


tgt 


highest peaks of the PHD, 


where [_■] denotes rounding to the nearest integer. 
Because of these two properties, the number of targets as well 
as their states can be estimated independently at each time 
step without any knowledge of their identities. That way, the 
data association issue is avoided. However, this also means 
that PHD Filters cannot deliver the continuous track of a 
specific target. If continuous tracks are required, an additional 
association step has to be performed. Two possible association 
algorithms for track continuity can be found in ED- 


D. The PHD Filter 

The PHD Filter is an approach for recursively propagating 
the PHD Djij(s(*)|Eo:i) at time step i given measurements up 
to time step i over time. If the RFS E is Poisson-distributed, 
then its PHD is equal to its intensity function and is, hence, 
a sufficient statistic 02) • In this case, the PHD recursion is 
given by the following prediction and update equations 02: 

A|i-l(s(*)l S 0:»-l) = &»(s(i)) + jPs{s{i - 1)) X 

fi\i-i(s(i)\s(i - l))A-i|i-i(«(* - l)|E 0: *_i)ds(i - 1) 

( 12 ) 


A|i(s(i)|£ 0: ») 


1 -PD + 


E 

zeSi 


_ Ppfi{z\s{i)) 

Afacfa(z) + p D J fi(z\s(i))D i \ i _ 


X -A|z— l) 


(a(*)|E 0: i_i)da(i) 

(13) 


Note that bi(s(i )) is the PHD of the birth set Bi of new 
targets appearing at time step i. In addition, ps(s(i — 1)) 
denotes the probability that a target survives the transition 
from time step i — 1 to i. The probability of survival de¬ 
pends on the previous state s(i — 1) because a target that 
is close to the border of the ROI and has a velocity vector 
pointing away from it is unlikely to be present at time step 
i. Furthermore, — 1)) and fi(z\s(i)) denote 

the transition probability and the likelihood, respectively. The 
probability of detection pd is constant since it is assumed that 
all targets can be detected if the ROI is covered. The term 
AfaCfa(z) represents Poisson-distributed false-alarms due to 
clutter, where Afa is the false alarm parameter, which is 
distributed according to its probability density function Cfa(z). 


III. Distributed Multi-Target Target Tracking 

In this section we introduce the Diffusion Particle PHD 
Filter (D-PPHDF), a distributed Particle Filter implementation 
of the PHD Filter for performing MTT in a sensor network 
without a fusion center. Before diving into the algorithm, we 
briefly review the concept of Adaptive Target Birth (ATB) and 
discuss the modification we applied in the D-PPHDF. 
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A. Adaptive Target Birth (ATB) 

Standard formulations of the PHD Filter consider the PHD 
bi(s(i)) of the birth set B, to be known a priori |f30l . For 
typical tracking applications such as air surveillance, this is a 
reasonable assumption since new targets should appear at the 
border of the ROI given continuous observation. An alternative 
is to make the target birth process adaptive and measurement- 
driven as suggested in OOl . Ufl . To this end, the PHD - and 
consequently the set of particles and weights approximating it 
in a Particle Filter implementation - is split into two densities 
corresponding to persistent objects, which have survived the 
transition from time step i — 1 to i, and newborn objects, 
respectively. 

In ED, OTI . the PHD of newborn objects is approximated 
by randomly placing Np new particles around each target 
measurement, with Np denoting the number of particles per 
target. We improve upon this approach by only considering 
measurements with no noticeable impact on any persistent 
particle weight, as these may indicate the appearance of a 
new target. That way, the number of newborn particles is 
further reduced and a possible overlap between persistent and 
newborn PHD is avoided. With the transition to time step 
i + 1, the newborn particles become persistent. Furtherore, 
we perform the ATB step towards the end of each iteration of 
the algorithm and only consider the particles representing the 
persistent PHD in the prediction, weighting, and resampling 
steps. Hence, the update equation (fl3t does not have to be 

modified as in ED, ED- 

While ATB delays the tracking algorithm by one time step, 
it is much more efficient as it only places new particles in 
regions in which a target is likely to be found. In addition, 
there is no need for an explicit initialization step since the 
first incoming target will trigger the deployment of a newborn 
particle cloud around its corresponding measurement. 

B. The Diffusion Particle PHD Filter (D-PPHDF) 

The proposed Diffusion Particle PHD Filter (D-PPHDF) is 
an extension of the single-sensor Particle PHD Filter (PPHDF) 
ED, ESI, Ea for the multi-sensor case. Furthermore, it 
relies on ATB for a more efficient target detection. The 
communication scheme we employ to exchange measurements 
and estimates between nodes is inspired by the two-step com¬ 
munication used in the context of Diffusion Adaptation 133]. 
However, the algorithm does not rely on least-mean-squares 
or any other kind of adaptive filter. First, each node k in the 
active part of the network obtains an intermediate estimate of 
the states of the targets present, i.e., of the PHD of persistent 

( 4 -^fc,pers(*) 

targets - represented by the set |Sfc jpers (*), <, pers (*) j 

of persistent particles with corresponding weights - based on 
neighborhood measurements. In other words, every active node 
runs a separate PPHDF with access to measurements from its 
neighborhood J\f k , defined as 

N k = {le{l,...,N}\ \\Xl - X k \\ 2 < Ream}, k=l,...,N., 

(14) 

where R com denotes the communication radius. Second, each 
active node combines the intermediate estimates from its 


neighborhood to a final, collaborative estimate. To this end, 
the persistent particle sets of all neighbors are merged into a 

collective set | sJ co11 (j),m;^ co11 (z) | of persistent neigh¬ 

borhood particles and corresponding weights before the clus¬ 
tering step, with Nk t coii(i) denoting the number of collective 
persistent neighborhood particles. In the sequel, we will look 
at the individual steps of the D-PPHDF in more detail: 

( 'j Afe )C0 it(f 1) 

. Merging: The sets i s p k n (i - 1), < u (i - 1) \ 

i ) p= i 

f 'l ttk, new(® 1) 

and (<new(* - 1 )»wjb, ne w(* ~ !)} =1 COnsist 0f 

the collective persistent neighborhood particles and 
newborn particles of node k, Sfc Coll (* — 1) and 
si (4 ^ 1 ), respectively, at time step i — 1 with their 
respective weights w p kco[x {i - 1) and vj p m Ji - 1). 
These sets are merged to become the total set 

{ s fc,tot(*).^fc,tot(*)/ of particles and weights 

of node k at time step i. Here, N klo t(i) is the total 
number of particles of node k at time step i, which is 
given by 

IVfc.tot(*) = N k , coll (i - 1) + lVfc,new(* - 1), (15) 

with iVfc con(z — 1) and N k , new (i — 1) denoting the re¬ 
spective number of persistent neighborhood and newborn 
particles at the previous time step. Note that since the 
sets of particles and weights represent PHDs, merging 
the sets corresponds to summing these PHDs. 

• Predicting: Each particle is propagated through the sys¬ 
tem model to become a persistent particle. The system 
model is assumed to be the same for each target and given 
by Equation <|TJ. Since the process noise is captured by 
the spread of the particle cloud, the respective term can 
be removed from the equation, yielding 

S fc,pers(*) = ^KtotW. P= !>■•■> (16) 

The corresponding weights are multiplied with the prob¬ 
ability of survival ps, which is assumed to be constant 
for the sake of simplicity according to 

<pers(^ - 1) = PSU% M (i), (1?) 

P — 1? • • • ? -^pers(^)* 

The prediction of particles and weights corresponds to 
the second term in Equation (fl2l) . 

• Measuring & Broadcasting (1): The sensor nodes obtain 
measurements of the targets and forward them to their 
neighbors. 

« Weighting: The persistent particle weights of node k are 
updated by applying a weighting step corresponding to 
Equation (IT3l i iteratively for each neighbor. Using the 
product operator, this weighting step can be compactly 
denoted as 

*A constant probability of survival ps is a reasonable assumption if the 
targets move relatively slowly with respect to the observation time and the 
size of the ROI. 
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l.persW = II 


1 ~PD + 


W k,j, update (*) 

Z) -es;. 


(18) 


w fc,per s (*K- !). 


with 


w 


s (0 = 


PDfi(Zj | a; p (i)) 


^OiUpdate\ ^ AfaCf A («j)+^(*j) : 


(19) 


where E* is the set of measurements obtained by node l 
and C{zj) is calculated as 

*Z lP _W 

£{zj)= PDMzj\x q (i))w q kpets {i\i-l). (20) 

9=1 


Note that ffzj \ x p (i)) is the likelihood and x p (i) is 
the location vector of particle p. 

Afterwards, each node k obtains the set E^ cand of can¬ 
didate measurements, i.e., measurements that are not 
responsible for the highest weighting of any persistent 
particle, to be used in the ATB step later on. The set 
E-'cand is found according to 


y \k 

cand 



m p = argmax w p kj} update (i), 

P = 1; •••) Nk, pers(t) j" ■ 

( 21 ) 


• Resampling: Each node k calculates its own expected 
number of targets N k ,igi{i) from its total persistent parti¬ 
cle mass according to 


Vfc,tgt(^) 


fc ,tot(t) 


H W fc,pers(*) 
P= 1 


( 22 ) 


Consequently, the number of persistent particles of node 
k is updated as 


N k 

,pers ) = N kM (i)N P . (23) 

Furthermore, the set of persistent particles of node k 
has to be resampled by drawing N kveis (i) particles with 
replacement from it. Note that the probability of drawing 
particle p is given by E ,1>eis since the weights do not 

-Nfc,tgt(*) 

sum to unity. Then, the weights are reset to equal values 
as 


^pers W 


N k , tgt (i) 

Nk, pers(^) 


P = lj.-.jJV^persW- ( 24 ) 


• Broadcasting (2): Every node k transmits its 
set of resampled persistent particles and weights 

( 'l IVfc,pers(t) 

( S fc,pers(*)»^fe,per s (*)j =1 to its neighbors. 

• Clustering: Each node k forms a collective set of persis¬ 
tent neighborhood particles s^. coll (i) and corresponding 


weights w p m (i) according to 


( s tV fe co u(i) 

{ s fc,cou(*)>^fc,cou(*)} p=1 


• if 'J N i t pers (i) ( 2 ^) 

U { S f,per s (*).<pers(*)} =1 > 


with 




N k , con (i) — ^ ' Ni tPers (i) 

l&Nk 


(26) 


denoting the number of collective persistent neighbor¬ 
hood particles of node k. As in the merging step, this 
corresponds to summing the corresponding PHDs to 
obtain an updated single-sensor PHD with a probabil¬ 
ity distribution reflecting the information of the entire 
neighborhood of node k. Note that the PHDs might not 
be independent if a target is detected by more than one 
neighbor. However, this is not a problem since merging 
the particle sets simply results in the respective target 
being represented by more particles. Hence, node k will 
be able to estimate the corresponding location more 
accurately. 

The estimated target states are found by clustering the 
collective persistent particles. Since the expected number 
of targets Ni jtgt (i), l £ J\f k might be different for each 
neighbor, we resort to hierarchical clustering of the 
single-linkage type |34| . Here, the sum of the expected 
number of targets over the neighborhood can serve as an 
upper bound for the number of clusters. Note, however, 
that if two targets are close to each other, clustering 
algorithms might not be able to resolve both targets 
correctly. 

• Roughening: A roughening step is performed to counter 
sample impoverishment [35.]. To this end, an independent 
jitter s 3 (i) is added to every resampled particle. Each 
component s 3 c (i), c = 1 of the jitter with di¬ 

mensionality d is sampled from the Gaussian distribution 
Af(0, (a 3 c (i)) 2 ). The component-wise standard deviation 
of the jitter is given by 

ai(i) = KE c N k , con (i)~ 1 / d , (27) 


where E c is the interval length between the maximum 
and minimum samples of the respective component. To 
avoid evaluating E c separately for each particle cluster, 
it is assigned an empirically found constant valued Note 
that d = 4 since the dimensionality of the jitter vector 
s 3 (i) and the particle state vector s p (i) have to coincide. 
In addition, K is a tuning constant, which controls the 
spread of the particle cloud. 

« Adaptive Target Birth: Np new particles are placed ran¬ 
domly around each candidate measurement Zj £ E^ cand 
leading to a total number of N k}tievl (i) = N P ■ |E^ cand | 
newborn particles for node k. Every newborn particle is 


2 Since the noise variances as well as the network topology are fixed, the 
true value of E c will not change significantly over time and between clusters, 
so this is a valid simplification. 










SUBMITTED FOR PUBLICATION 


6 


associated with a weight that is chosen according to 


w 


P 

/c,new 


(<) = 


Pb 


^Lew«’ 


P= l,...,A r fc,new(«), (28) 


where p/j is the probability of birth. Depending on the 
application, pp can depend on time as well as on the 
location of the respective particle. For simplicity, the 
probability that a new target enters the ROI is assumed 
to be equal for all locations in the birth region over time. 
The target birth process corresponds to the first term in 
Equation (IT2ll . 

Figure ITbl ) shows an example of tracking three targets, 
which move along the deterministic tracks depicted in Fig¬ 
ure flat . using the D-PPHDF. Note that each small colored 
dot corresponds to a target location estimate obtained by the 
respective node with the same color while the light grey dots 
represent the collective measurements from all nodes. From 
this illustration, the following properties of the D-PPHDF are 
apparent: First, the algorithm only delivers separate location 
estimates - represented by the small colored dots - for 
each time instant rather than continuous tracks, which - as 
mentioned before - is a common property of PHD filters. 
Second, the network as a whole would be able to correctly 
track all three targets, while a single node only obtains the 
locally relevant subtracks of the targets in its vicinity. Third, 
the employed two-step communication theme is able to extend 
the vicinity of a node far beyond its own sensing radius of 
R sen = 6 m. This can, for instance, be seen from the fact that 
the lime-green node located at [-14,-23] is able to obtain 
location estimates of target 2, which enters the ROI from 
the south. Finally, Figure [Tbl also illustrates the resolution 
problem of clustering. When targets 1 and 2, which enter the 
ROI from the north and the south, respectively, cross paths, the 
nodes in their vicinity see them as just one target. This leads 
to an aggregation of target location estimates around [9, 0]. 

The pseudo-code of the D-PPHDF is given in Table UU 


C. Computational Complexity and Communication Load 

In this section we take a look at the computational com¬ 
plexity and the communication load the D-PPHDF imposes 
on each node in the active subnetwork. The following steps 
are performed at every time instant i but time dependency is 
omitted for simplicity. Note that each of the steps scales with 
the number of active nodes when considering the computa¬ 
tional complexity of the network as a whole. 

• Prediction: The prediction step described by Equations 
( IT6l l and (ITTIl is performed for each particle at every active 
node. Hence, it scales with the number of particles Nk,tot 
and the dimensionality d of the particle vectors. In order 
to obtain a tractable expression for the computational 
complexity, we assume each node to have the same 
number of particles N tot . 

=► 0(N tot d) 

• Weighting: Each particle is updated in the weighting step 
given by Equations ( fl8l - d2ll . The weight update as well 
as the designation of candidate measurements for ATB 
depends on the neighborhood size | A41 of node k and the 


number of measurements |E ; | of each of its neighbors l. 
For tractability reasons, we assume each node to have the 
same number of neighbors 2V n b and to obtain the same 
number of measurements iV me as- 
=F 0(N tot N ab N meils ) 

• Resampling: The estimation of the number of targets and 
the resampling step in Equations d22l) - (l24l) are linear in 
the number of particles used for the calculation |36|. For 
the sake of simplicity, we assume each active node to 
have the same estimate of the number of targets IV t g t . 

0(N t ot + WactiyeiVtgtiVp) 

• Clustering: The complexity of single-linkage clustering 
is cubic in the number of particles, i.e., in the number 
of neighbors N nb of each node, the estimated number of 
targets N lgl , the number of particles per target Np, and 
the dimensionality d of the particles lf37l . 

=► 0({N nb N tgt N P d) 3 ) 

m Roughening: Roughening (Equation 1271 )1 is performed 
for every collective particle and is linear in the dimen¬ 
sionality of the particles. 

=► 0(N nb N tgt N P d) 

• Adaptive Target Birth: The birth process depends on the 
number of particles per target Np as well as the number 
of candidate measurements A ? can( i, which is assumed equal 
for each active node to ensure tractibility. 

=F 0(N P N can d ) 

As far as the communication load is concerned, the D-PPHDF 
requires the broadcasting of measurements, i.e., 2 scalars 
per measurement, over the neighborhood in the first broad¬ 
casting step. In the second step, the sets of particles and 
weights, i.e., 5 scalars per particle, are transmitted. Clearly, 
the communication load strongly depends on the number of 
nodes in the network, or more precisely the number of active 
nodes and the size of their respective neighborhood. As an 
extension of the D-PPHDF, one could think of changing 
the second broadcasting step and transmit Gaussian Mixture 
Model representations - instead of the actual particles and 
weights - that will be resampled at the receiver node. That 
way, communication load could be reduced to transmitting 
only a few scalars in the second broadcasting step at the cost of 
estimation accuracy and additional computational complexity. 
However, a thorough treatment of this extension is beyond the 
scope of this work. 

IV. Centralized Multi-Target Tracking 

Having presented the D-PPHDF as a distributed solution 
for MTT in a sensor network, we propose the centralized 
counterpart to our approach in the sequel. 

A. The Multi-Sensor Particle PHD Filter (MS-PPHDF) 

The proposed Multi-Sensor Particle PHD Filter 
(MS-PPHDF) is a centralized, multi-sensor PPHDF that 
relies on a fusion center with access to the measurements 
of all nodes in the network. It is based on the formulation 
of the single-sensor PPHDF in ED, EH, 01 but with an 
extended measurement set comprising the measurements of 
the entire network. Hence, one might obtain more than one 
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measurements per target - a change to the typical assumption 
in target tracking that each target produces at most one 
measurement ED. To account for this change, we add a 
pre-clustering step before the weighting step and normalize 
the weight update accordingly. A similar partitioning of the 
measurement set is used in extended target tracking, where 
a sensor can receive multiple target reflections due to the 
target’s physical extent (38), (39) . 

In the following, we will look at the individual steps of the 
algorithm in more detail: 

. Merging-. The sets {s^ ra (z - 1), w$ eTS (i - l)}^ (l_1) 

and {snew(* — l),zt;ne W (z — 1) cons i st °f the 

persistent and newborn particles, Sp ers (« —1) and Sn ew (* — 
1), respectively, at time step i — 1 with their respective 
weights wipers(i — 1) and Wnew{i — 1)- These sets are 
merged to become the total set {stot(*)> w tot(*)}^5'i^ of 
particles and weights at time step i. Here, is the 

total number of particles at time step i, which is given 

by 


N tot (i) = N peTS (i - 1) + N nev/ (i - 1), (29) 

with iVpers (i — 1) and N new (i — 1) denoting the respec¬ 
tive number of persistent and newborn particles at the 
previous time step. 

• Predicting: As in the D-PPHDF, each particle is propa¬ 
gated through the system model according to 

S pers(*) = -Fs?ot(*), P= -Npers(i) = N tot (i) 

(30) 


to become a persistent particle. The corresponding 
weights are multiplied with the probability of survival 

PS as 

<ers(^ ~ 1) = PSV&tii), 

P — 1? • • • 5 -^pers(^)* 

• Measuring: The sensor nodes obtain measurements of the 
targets. 

• Pre-Clustering: Since there might be more than one 
measurement per target, the measurements of the entire 
network are pre-clustered before the weighting step and 
each measurement is assigned a label C(z) that reflects 
the cardinality of its own cluster. This can be done, for 
instance, using single-linkage clustering (34) . The clus¬ 
tering is based on the distance between measurements, 
i.e., spatially close measurements are assumed to stem 
from the same target. Hence, when two or more targets 
are too close to each other, cardinality errors may occur. 

• Weighting: All available target measurements, which 
comprise the set £*, are used to update the persistent 
particle weights according to 


«&»(*) = 


1 ~PD+ Y W l updated) 


ZjGSi 


<ers(*l* - 1), 


with 


p D fi(zj | xP(i)) 


w p (t) = _ 

,,update! J (A FA C FA (2 j) + C{Zj)) C(Zj) 


(32) 

, (33) 


and 


TVpers (») 

C(Zj) = Y PdM z j I ® 9 (*)Vpers(*l* - !)• (34) 

9=1 

Note that - in contrast to the D-PPHDF - the weighting 
step is applied only once using the entire set of mea¬ 
surements. Therefore - and since there might be more 
than one measurement per target - we have to ensure 
that the weight update terms w? update - and consequently 
the particle weights - still sum to the number of targets 
present. This is done by normalizing Equation (f33l) with 
C(zj), i.e., the cardinality of the cluster to which the 
current measurement Zj belongs. 

Afterwards, we form the set £i,cand of candidate measure¬ 
ments for the ATB step according to 


-'i.eand 


= E, 


m p = argmax w£ update (i), 

J (35) 


P — 1 1 ■ • ■) -^pers(*)}• 


• Resampling: The expected number of targets N tgt (i ) is 
calculated from the total persistent particle mass as 




W pers(*) 


P= 1 


(36) 


Consequently, the number of persistent particles is up¬ 
dated according to 


N peIS (i) =N tgt (i)N P . (37) 


Furthermore, the set of persistent particles is resampled 
by drawing N peis (i) particles with probability 


„(* 


Then, the weights are reset to equal values as 


IV, g ,(i) 


<ers00 = P=l,...,N pets {l). (38) 

m Clustering: In contrast to the D-PPHDF, there is only one 
estimate of the expected number of targets. Hence, we can 
use k-means clustering (40l to find the estimated target 
states by grouping the resampled particles into N tgt (i) 
clusters and calculating the centroid of each cluster. 

• Roughening: Roughening is performed analogously to the 
D-PPHDF. 

• Adaptive Target Birth: Np new particles are placed ran¬ 
domly around each candidate measurement z :] € £i jC and 
yielding a total number of N nev/ (i) = Np ■ |£i iC and| 
newborn particles. The corresponding weights are chosen 
according to 

W new(») = p = l,...,N asw (l), (39) 

where pp is the probability of birth. 

The pseudo-code of the MS-PPHDF is given in Table Hill 
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B. Computational Complexity and Communication Load 

In this section we analyze the computational complexity 
and the communication load of the MS-PPHDF. The following 
steps are performed at every time instant i but time dependency 
is omitted for simplicity: 

• Prediction: The prediction step described by Equations 
i l30ll and (OTb is performed for each of the iV to t particles 
and is linear in the dimensionality d. 

=> O(N t 0 t d) 

m Pre-Clustering: The pre-clustering step relies on single¬ 
linkage clustering. The complexity is therefore cubic in 
the total number of measurements N meas . E3 

=► O(^meas) 

• Weighting: Each particle is updated in the weighting step 
given by Equations (|321)-(|34||. The weight update as well 
as the designation of candidate measurements for ATB 
depends on the number of measurements N meas = |S|. 
=>■ ©(-/VtotiVmeas) 

• Resampling: The estimation of the number of targets and 
the resampling step in Equations l l36l) - (l38l) are linear in 
the number of particles used for the calculation If36ll . 

=> O{N t0t + N tgt N P ) 

m Clustering: In contrast to the D-PPHDF we can use 
fc-means clustering. The complexity of Lloyd’s imple¬ 
mentation is given by ED 
=> 0((N tgt Np) dN, » +1 \og(N tgt Np)). 

• Roughening: Roughening is linear in the dimensionality 
of the particles and their number. 

=> 0(N tgt N P d ) 

• Adaptive Target Birth: The birth process depends on the 
number of particles per target Np as well as the number 
of candidate measurements N cm d. 

=> (Di^Np iVcand) 

In summary, the computational complexity of the 
MS-PPHDF is largely comparable to that of the D-PPHDF. 
The only exception is the pre-custering step, which scales 
cubicly with the total number of measurements and adds 
additional complexity to the algorithm. As a tradeoff the com¬ 
munication load of the MS-PPHDF clearly is lower compared 
to the D-PPHDF because there is only the initial transmission 
of measurements from the nodes to the fusion center. However, 
considering a setup with relatively small communication radii, 
this initial communication step requires a lot of relaying and 
leads to high traffic density in the vicinity of the fusion center. 
Furthermore, this communication structure exhibits a single 
point of failure while a distributed sensor network is inherently 
redundant. 


V. Simulations 

In this section, we evaluate the performance of the 
D-PPHDF as well as the MS-PPHDF for tracking multiple 
targets in a sensor network with 1-coverage. In addition, we 
compare them to the alternative distributed PPHDF from lf20l . 
which will be referred to as Distributed Data Fusion Particle 
PHD Filter (DDF-PPHDF). Here, each node runs its own 
PPHDF using only its own measurements. In a subsequent 
step, the particles are distributed over the neighborhood and 


reweighted by fusing their corresponding Exponential Mixture 
Densities. 

Furthermore, we formulate the DPCRLB as a lower bound 
for evaluating the performance of the three algorithms in terms 
of the OSPA H5l metric. In our simulations, we compute 
the OSPA metric with respect to the joint set of target state 
estimates of the entire active network. The latter is found 
by clustering the target state estimates of all active nodes. 
Furthermore, we consider the squared OSPA metric scaled by 

the number of targets, i.e., jY tgt • ^ dp C ^ , as in lf42l . That way, 
we can use the DPCRLB, which will be introduced in the 
following, as a benchmark. 


A. The Distributed Posterior Cramer-Rao Lower Bound 
(DPCRLB) 


Rather than evaluating the performance of the different 
MTT algorithms based on an error metric, it makes more sense 
to derive a minimum variance bound on the estimation error, 
which enables an absolute performance evaluation. For time- 
invariant statistical models, the most commonly used bound 
is the Cramer-Rao Lower Bound (CRLB), which is given by 
the inverse of Fisher’s information matrix E2l. In El and 
|44| . the CRLB is used in the context of multi-sensor MTT of 
an unknown number of unlabeled targets in order to evaluate 
the performance, as well as prove the asymptotic efficiency 
of the PHD as the number of nodes goes to infinity. Since 
we are more interested in the tracking behavior of a fixed 
network over time, we resort to the Posterior Cramer-Rao 
Lower Bound (PCRLB), which is an extension of the CRLB 
for the time-variant case ESI. This bound can be calculated 
sequentially with the help of a Riccati-like recursion derived 
in El- Furthermore, in E71 and ESI , the PCRLB is adapted 
for an MTT scenario in which the tracker can obtain more 
than one measurement per target. 

Let 7 r™,m = 1 denote the probability that any 

measurement is associated with target m at time instant i 
as defined in 01- With the corresponding stochastic process 
II" 1 , the new stochastic process of association probabilities and 
target states to be estimated becomes = ( H} :M , S ( 1M J. 


Fisher’s information matrix J# i = 



Jnl 


can now be 


formed as described in E6l and fl7l . However, as the number 
of targets varies over time, i.e., targets might enter or exit the 
ROI, ./#, has to be expanded or shrunk in the inverse matrix 
domain as described in E71 . The PCRLB H, at time instant 
i can be obtained as the trace of the inverted submatrix ./=, 
according to |46j 



jW T-l 

J 3i J Yli 




(40) 


Note that, in a distributed MTT scenario, Bi corresponds to 
a lower bound on the estimation error of a central processing 
unit with access to all measurements. Since we are interested 
in completely distributed MTT with in-network processing, we 
extend the PCRLB to the Distributed Posterior Cramer-Rao 
Lower Bound (DPCRLB). To this end, each node k computes 
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its own PCRLB Bf considering only the measurements of its 
two-hop neighborhood, which is given by 

NV = U ( 4 D 

l&Mk 

i.e., the neighbors of node k and their neighbors. Furthermore, 
only the targets within the sensing range of A4 are taken 
into account. Clearly, only nodes with a neighborhood in 
the vicinity of at least one target will be able to calculate a 
PCRLB. The DPCRLB at time instant i is then obtained 
by averaging over these values according to 

U-) 

1 1 fcGVW 

where M is the set of all nodes that are able to compute a 
PCRLB. 

B. Simulation Setup 

In the following simulations, a static sensor network as 
depicted in Figure [Taj) is used to perform MTT. The network 
is centered around the point of origin [0, 0] T and distributed 
such that 1-coverage of the ROI is guaranteed. It covers an 
area of approximately 2500 m 2 . Clutter is assumed to be 
uniformly distributed over the sensing range of each node 
with an average rate of Afa = 0.1. For the sake of simplicity, 
collisions between targets and sensor nodes are neglected. An 
overview of all simulation parameters is given in Table UJ Since 
the purpose of this work is to introduce the MS-PPHDF as well 
as the D-PPHDF, verify their functionality, and compare them 
to alternative approaches, we consider a rather simple scenario 
with a high probability of detection and a low clutter level. In 
our future work, we will study more sophisticated scenarios 
to define possible breakdown points of our algorithms. 

We use the MS-PPHDF, the D-PPHDF, as well as the 
DDF-PPHDF to track three targets for i = 0,30. The 
targets enter the ROI at time steps * = 0,9, 14 from the north, 
south, and west, respectively. A Monte Carlo simulation with 
Nmc = 1000 runs is performed to evaluate the performance 
of the tracking algorithms in terms of the Optimal Subpattern 
Assignment (OSPA) metric and compute the DPCRLB as a 
benchmark. Note that the target trajectories as shown in Fig¬ 
ure [Ta)) are deterministic, as is often the case in target tracking 
simulations {471 in order to guarantee the comparability of the 
different Monte Carlo runs regarding, for instance, the number 
of targets present. 

C. Simulation Results 

The simulation results for tracking three targets in the given 
setup are depicted in Figure [2] While the top part considers 
zero-mean Gaussian measurement noise with a per-component 
variance of <r 2 = 0.1, the bottom part shows the results 
for cr 2 = 0.3. In addition to evaluating the performance of 
MS-PPHDF, D-PPHDF, and DDF-PPHDF in terms of the 
squared and scaled OSPA metric over time and comparing it 
to the DPCRLB as can be seen in the left part of the figure, we 
also look at the estimated number of targets, which is depicted 


TABLE I: Simulation parameters 


Variable 

Value 

Description 

Az 

1 S 

time step of the tracking algorithm 

N 

30 

number of nodes 

al 

0.1, 0.3 m 2 

componentwise power of meas. noise 

°\ 

0.01 m 2 

componentwise power of state noise 

-^com 

2i?sen 

communication radius 

Rsen 

6 m 

sensing radius 

E c 

6 

empirical interval length for jitter 

K 

0.2 

tuning constant for roughening 

N P 

500 

number of particles per target 

Pb 

0.8 

probability of birth 

Pd 

0.95 

probability of detection 

PS 

0.98 

probability of survival 

^FA 

0.1 

average no. of false alarms / clutter 

cfa(z) 

1 

TT-Rsen 

PDF of false alarms / clutter (uniform) 

c 

2 

cut-off value (OSPA) 

p 

2 

order of the OSPA metric 


in the right column. Since the OSPA metric contains a penalty 
for an erroneous estimate of the number of targets, this side- 
by-side comparison facilitates the interpretation of the tracking 
results. 

Let us start by considering Figures [20 and 1 2 hi ), i.e., the 
case of cr 2 = 0.1. First of all, we observe that neither 
tracking algorithm provides an OSPA value or an estimate 
of the number of targets for i = 0. This is expected and 
due to ATB, which initializes new particle clouds based on 
the measurements from the previous time step. Thus, target 
birth is delayed by one time step and tracking can only be 
performed for i > 0. The same effect can be witnessed at 
i = 9 and i = 14, respectively, which mark the time instants 
at which targets 2 and 3 enter the ROI. Here, the OSPA curves 
of all trackers exhibit a spike, which is due to the fact that the 
newborn particles are not yet considered in the tracker and, 
hence, the number of estimated targets is too low, as can be 
seen in Figure [201. 

Another sudden rise of all the OSPA curves can be observed 
in the time interval 20 < i < 24 with a valley at i = 22. 
Looking at the estimated number of targets, we can attribute 
this phenomenon to the fact that only two of the three targets 
are recognized by the tracking algorithms. Since the target 
trajectories are deterministic, we know that in the given time 
interval targets 2 and 3 cross paths. Due to the inability of 
the clustering algorithm to separate strongly overlapping sets 
of measurements, the two targets merge into one as long as 
they are close to each other. When the two targets occupy 
almost exactly the same position, i.e., at i = 22, the OSPA 
metric decreases due to the decrease in measurement variance. 
As the targets drift apart, the variance and with it the OSPA 
metric increases up to the point where the two targets can be 
recognized as separate again and the corresponding penalty is 
switched off. 

Looking at the overall picture in Figure [20, which shows 
the case of cr 2 = 0.1, it is evident that the centralized 
MS-PPHDF and the distributed D-PPHDF achieve approx¬ 
imately the same performance with OSPA values closely 
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(a) Squared and scaled OSPA compared to DPCRLB, rC =0.1 




(c) Squared and scaled OSPA compared to DPCRLB, cr^ = 0.3 



Fig. 2: Simulation results for tracking three targets in the ROI over time using the Multi-Sensor Particle PHD Fil¬ 
ter (MS-PPHDF), the Diffusion Particle PHD Filter (D-PPHDF), and the Distributed Data Fusion Particle PHD Filter 
(DDF-PPHDF). The left part of the figure shows the squared and scaled Optimal Subpattern Assignment (OSPA) metric 
for each algorithm compared to the Distributed Posterior Cramer-Rao Lower Bound (DPCRLB), while the right part compares 
the estimated to the true number of targets. 


approaching the DPCRLB when the number of targets stays 
constant. Furthermore, both algorithms deliver very accurate 
estimates of the number of targets, given they are separable 
by clustering, as can be seen in Figure l2bt. The DDF-PPHDF, 
however, continuously exhibits a worse performance than the 
D-PPHDF, both in terms of the OSPA metric as well as 
the estimated number of targets. This is where the additional 
communication in the proposed D-PPHDF shows its strength 
in reducing uncertainty due to measurement noise and clutter. 
Apart from achieving worse tracking results, the DDF-PPHDF 
also has more difficulty in separating targets 1 and 2 when they 
cross paths, resulting in an earlier rise and a later fall of the 
OSPA metric, compared to our approach. 

In the case of of = 0.3, the overall performance of the 
different tracking algorithms is very similar to the case of 
of = 0.1. In order to make a statement on how the different 


tracking algorithms compare, let us neglect the penalty due 
to an erroneous estimate of the number of targets and take a 
look at Figures [3al> andl3bl). which are zoomed-in versions of 
Figures l2all and [2c1 ). respectively. 

In Figures[3a]i andITbli the DPCRLB is given as a benchmark 
for tracking performance. One can observe that its value 
is always smaller or equal to the respective measurement 
variance. As stated before, the centralized MS-PPHDF and 
the distributed D-PPHDF exhibit very similar performance and 
deliver better tracking results than the DDF-PPHDF. While the 
MS-PPHDF achieves lower OSPA values than the D-PPHDF 
when the number of targets stays constant, i.e., for 3 < i < 8 
and 24 < i < 30, the D-PPHDF performs better directly after 
a new target appears, i.e., for 1 < i < 2, 10 < * < 13, and 
15 < i < 18. This is likely due to the fact that the two-step 
communication scheme employed in the D-PPHDF is able to 
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(a) Squared and scaled OSPA compared to DPCRLB, cr^ =0.1 (b) Squared and scaled OSPA compared to DPCRLB , = 0.3 

Fig. 3: Simulation results for tracking three targets in the ROI over time (zoomed in). The squared and scaled Optimal Subpattern 
Assignment (OSPA) metrics of using the Multi-Sensor Particle PHD Filter (MS-PPHDF), the Diffusion Particle PHD Filter 
(D-PPHDF), and the Distributed Data Fusion Particle PHD Filter (DDF-PPHDF), are compared to the Distributed Posterior 
Cramer-Rao Lower Bound (DPCRLB). 


reduce the impact of measurement noise and clutter faster than 
the centralized MS-PPHDF can. 

Looking at the case of crj =0.3 in Figure [3bt. we observe 
that the higher measurement noise affects the performance 
of all algorithms, resulting in higher OSPA curves. While 
the OSPA curves of the MS-PPHDF and the DDF-PPHDF 
are proportionally shifted upward by approximately the same 
value, i.e., they are equally impacted by the higher noise 
level, the D-PPHDF seems to be slightly more affected by 
the change. However it still outperforms the DDF-PPHDF at 
all time instants. 

All in all, the proposed D-PPHDF yields better performance 
than the existing DDF-PPHDF in estimating the number of 
targets and tracking them, irrespective of the amount of mea¬ 
surement noise. In addition, it is also a bit faster in delivering 
correct state estimates of new targets than the centralized 
MS-PPHDF and performs only slightly worse once the number 
of targets stays constant. In our future work, we will look at 
ways to further improve the performance of the MS-PPHDF 
and the D-PPHDF in order to approach the DPCRLB even 
more closely. 

VI. Conclusion 

In this work, we developed a distributed as well as a 
centralized PPHDF for MTT in sensor network. We, further¬ 
more, came up with a distributed version of the PCRLB that 
served as a benchmark in the performance evaluation. Our 
simulation results showed that the distributed D-PPHDF is 
faster in correctly tracking new targets than the centralized 
MS-PPHDF and performs only slightly worse when the num¬ 
ber of targets stays constant. In addition, it delivers accurate 
tracking results as long as the targets are far enough apart 
so that their corresponding measurement clouds are separable. 
Our approach outperforms the existing DDF-PPHDF at the 
cost of additional communication between sensor nodes. 
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1: input: d, E c , K, n, N, Np,pp, p$, Afa> cfa 

2: initialize: {^, coll (0), «>^ coll (0)}^’ coll(0) = |^ new (0), ^ new (0)}^’ new(0) =0. 

3: while i < n do 
4: for k = 1,..., N do 

5: Merge the sets of collective persistent and newborn particles with corresponding weights: 


(s p ,(0) ’ =(s? n (i-l)V ' coll( J u{< (i-1)} *■' 

t k, tot'- n k, tot'- ' J p—i l /c,coll v k, colD ' J p = \ l k, new' n k, new' ' J p—1 


N k, new!*- 1 ) 


6: for p = 1,..., N kM (i) do 

7: Predict the new state of each particle and update the weight with the probability of survival pg: 

S fc,persM = Fs fc,totW 


" ; fc,per S ( i H- 1 )=PS<, tot »- 


Update the weights using neighborhood measurements: 


sW= II 1 -PD+ update^ <,per S (*l* “ *)> 

1 C- Kf . - _ z- v>7 


W l. Mandate (*) = 


PDfi(zj I * P W) 


fc,i,update^ Afacfa(z . ) + £(z .)’ 

N k,fas(i) 

£( z i)=* ^2 PDfi{zj\x q (i))wl ptts (i\i-l). 

9=1 


9: end for 

10: Form the set of candidate measurements for ATB: 


S £cand = S ?\{ z m p m P = ar § “’fcj,update W> P = 1, JVfc.persW, j = L •••> |^ | VZ € Affc }. 


11: Calculate the estimated number of targets: 


A,tg .» = X! ^,pers(*) • 

P=1 

12: Resample iV/ C)Pers (i) = Nk,igii})Np particles and reset the weights: 

w p (i) = 1 N (i) 

fc ’ peK JVfc.persW’ ,PeIS ' ' 

13: Merge the sets of persistent neighborhood particles and weights: 

{<,con(*)t^, co u(*)}"l fc r (i) = U { S r,per s (*)t<pe rs (*)}!lr (i) 


14: Use single-linkage clustering to identify iVtgt(i) clusters and find the set of estimated target states by calculating the centroids 

15: Add an independent jitter to each particle using a component-wise standard deviation of: 

ai(i) = KE c N ktCO \i(i)— 1 t d . 

16: Place Np new particles randomly around each candidate measurement Zj E cand . Set the weights as: 


v r\ rts 

W, (2 ) = -, 

fc,new v 1 at /-\ ’ 


P — ^fc,new(*)‘ 


end for 

i <— i + 1 

end while 
return 


TABLE II: The Diffusion Particle PHD Filter. 
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1: input: d, E c , K, n, N, Np,pp, ps, Afa, cfa 

2: initialize: {s^ oU (0), tnf oll (0)} p ”“ (0) = {s£ ew (0), «’new(0)}^“* (0) =0. 

3: while i < n do 

4: Merge the sets of persistent and newborn particles with corresponding weights: 


= {« - l),<=rs(* - ' U {^(i ~ l),<w(* - 1)}^ 


Nnew (* 1) 


5: for p = 1, ..., JVtot(i) do 

6: Predict the new state of each particle and update the weight with the probability of survival ps : 


«pers(*) = Fs? ot (*) 
«’pers(*l*- 1) =PSW\ m(t). 


7: Cluster the measurements using single-linkage clustering and assign each measurement 2 : a label C(z) reflecting the cardinality of its cluster. 

8: Update the weights using the measurements of the entire network: 


^per S (*) = 1 ~PD+ ^2 updated) “ U 

^ .CV ■ 


w v (i) = PdMzj I * p (0) 

i.upda'eW (Afacfa(z . )+£(z . ))C(z . ) ’ 

Npers(i) 

£(zj) = ^2 Pofi( z j I x q (i))w$ as (i\i - 1). 
q =1 


9: end for 

10: Form the set of candidate measurements for ATB: 


Si.cand = £j\|z mp rn p = argmax w P update («), P = 1, N fas (i), 3 = 1, •••, |£*| j. 


11: Calculate the estimated number of targets: 


AW*) 

Agt(*) = ^2 TO &ra(®) • 

_ p= 1 


12: Resample 7Vp er s(i) = Ntgt(i)Np particles and reset the weights: 


<rs« = P = 1, • • • , ^pers(i). 


Apers (i) 


13: Find the set of estimated target states by using fc-means clustering and calculating the centroid of each cluster. 

14: Add an independent jitter to each particle using a component-wise standard deviation of: 

a{(i) = K E c N fe[s (i)— 1 / d . 

15: Place Np new particles randomly around each candidate measurement Zj G S^cand- Set the weights as: 


tttnew(^) 


Wnew(*) ’ 


P — -^new(^)* 


i <— i + 1 

end while 
return 


TABLE III: The Multi-Sensor Particle PHD Filter. 
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